// complex standard header -*-c++-*-
// Copyright 2003-2010 IAR Systems AB. 
#ifndef _COMPLEX_
#define _COMPLEX_

#ifndef _SYSTEM_BUILD
#pragma system_include
#endif

#include <ymath.h>
#include <ccomplex>
#include <cmath>
#include <sstream>

_C_STD_BEGIN
 #ifndef _C_COMPLEX_T
  #define _C_COMPLEX_T

typedef struct _C_double_complex
        {       /* double complex */
        double _Val[2];
        } _C_double_complex;

typedef struct _C_float_complex
        {       /* float complex */
        float _Val[2];
        } _C_float_complex;

typedef struct _C_ldouble_complex
        {       /* long double complex */
        long double _Val[2];
        } _C_ldouble_complex;
 #endif /* _C_COMPLEX_T */

_C_STD_END

        // COMPLEX _Val OFFSETS
 #define _RE    0
 #define _IM    1

_STD_BEGIN
typedef _CSTD _C_double_complex _Dcomplex_value;
typedef _CSTD _C_float_complex _Fcomplex_value;
typedef _CSTD _C_ldouble_complex _Lcomplex_value;

#define __STD_COMPLEX   /* signal presence of complex classes */


                // CLASS _Ctraits_double
class _Ctraits_double
{       // complex traits for double
public:
  typedef double _Ty;

  static _Ty _Cosh(_Ty _Left, _Ty _Right)
  {       // return cosh(_Left) * _Right
    return (_CSTD cosh(_Left) * _Right);
  }

  static short _Exp(_Ty *_Pleft, _Ty _Right, short _Exponent)
  {       // compute exp(*_Pleft) * _Right * 2 ^ _Exponent
    return (_CSTD exp(*_Pleft) * _Right * (_Ty)(1 << _Exponent));
  }

  static _Ty _Infv(_Ty)
  {       // return infinity
    return 0.Infinity;
  }

  static bool _Isinf(_Ty _Left)
  {       // test for infinity
    return isinf(_Left);
  }

  static bool _Isnan(_Ty _Left)
  {       // test for NaN
    return isnan(_Left);
  }

  static _Ty _Nanv(_Ty)
  {       // return NaN
    return 0.NaN;
  }

  static _Ty _Sinh(_Ty _Left, _Ty _Right)
  {       // return sinh(_Left) * _Right
    return (_CSTD sinh(_Left) * _Right);
  }

  static _Ty atan2(_Ty _Yval, _Ty _Xval)
  {       // return atan(_Yval / _Xval)
    return (_CSTD atan2(_Yval, _Xval));
  }

  static _Ty cos(_Ty _Left)
  {       // return cos(_Left)
    return (_CSTD cos(_Left));
  }

  static _Ty exp(_Ty _Left)
  {       // return exp(_Left)
    return (_CSTD exp(_Left));
  }

  static _Ty ldexp(_Ty _Left, int _Exponent)
  {       // return _Left * 2 ^ _Exponent
    return (_CSTD ldexp(_Left, _Exponent));
  }

  static _Ty log(_Ty _Left)
  {       // return log(_Left)
    return (_CSTD log(_Left));
  }

  static _Ty pow(_Ty _Left, _Ty _Right)
  {       // return _Left ^ _Right
    return (_CSTD pow(_Left, _Right));
  }

  static _Ty sin(_Ty _Left)
  {       // return sin(_Left)
    return (_CSTD sin(_Left));
  }

  static _Ty sqrt(_Ty _Left)
  {       // return sqrt(_Left)
    return (_CSTD sqrt(_Left));
  }

  static _Ty tan(_Ty _Left)
  {       // return tan(_Left)
    return (_CSTD tan(_Left));
  }
};

// CLASS _Ctraits_float
class _Ctraits_float
{       // complex traits for float
public:
  typedef float _Ty;

  static _Ty _Cosh(_Ty _Left, _Ty _Right)
  {       // return cosh(_Left) * _Right
    return (_CSTD _F_FUN(cosh)(_Left) * _Right);
  }

  static short _Exp(_Ty *_Pleft, _Ty _Right, short _Exponent)
  {       // compute exp(*_Pleft) * _Right * 2 ^ _Exponent
    return (_CSTD _F_FUN(exp)(*_Pleft) * _Right * (_Ty)(1 << _Exponent));
  }

  static _Ty _Infv(_Ty)
  {       // return infinity
    return 0.Infinity;
  }

  static bool _Isinf(_Ty _Left)
  {       // test for infinity
    return isinf(_Left);
  }

  static bool _Isnan(_Ty _Left)
  {       // test for NaN
    return isnan(_Left);
  }

  static _Ty _Nanv(_Ty)
  {       // return NaN
    return 0.NaN;
  }

  static _Ty _Sinh(_Ty _Left, _Ty _Right)
  {       // return sinh(_Left) * _Right
    return (_CSTD sinhf(_Left) * _Right);
  }

  static _Ty atan2(_Ty _Yval, _Ty _Xval)
  {       // return atan(_Yval / _Xval)
    return (_CSTD _F_FUN(atan2)(_Yval, _Xval));
  }

  static _Ty cos(_Ty _Left)
  {       // return cos(_Left)
    return (_CSTD _F_FUN(cos)(_Left));
  }

  static _Ty exp(_Ty _Left)
  {       // return exp(_Left)
    return (_CSTD _F_FUN(exp)(_Left));
  }

  static _Ty ldexp(_Ty _Left, int _Exponent)
  {       // return _Left * 2 ^ _Exponent
    return (_CSTD _F_FUN(ldexp)(_Left, _Exponent));
  }

  static _Ty log(_Ty _Left)
  {       // return log(_Left)
    return (_CSTD _F_FUN(log)(_Left));
  }

  static _Ty pow(_Ty _Left, _Ty _Right)
  {       // return _Left ^ _Right
    return (_CSTD _F_FUN(pow)(_Left, _Right));
  }

  static _Ty sin(_Ty _Left)
  {       // return sin(_Left)
    return (_CSTD _F_FUN(sin)(_Left));
  }

  static _Ty sqrt(_Ty _Left)
  {       // return sqrt(_Left)
    return (_CSTD _F_FUN(sqrt)(_Left));
  }

  static _Ty tan(_Ty _Left)
  {       // return tan(_Left)
    return (_CSTD _F_FUN(tan)(_Left));
  }
};

class double_complex;

// CLASS float_complex
class float_complex
  : public _Fcomplex_value
{       // complex with float components
public:
  typedef float _Ty;
  typedef float_complex _Myt;
  typedef _Ctraits_float _Myctraits;
  typedef _Ty value_type;
  // Used to get xcomplex 'templates' correct.
  typedef _Fcomplex_value _MyBaseT;

  explicit float_complex(const double_complex&);  // defined below

  float_complex(const _Ty& _Realval = 0, const _Ty& _Imagval = 0)
  {       // construct from float components
    _Val[_RE] = _Realval;
    _Val[_IM] = _Imagval;
  }

  float_complex(const _Fcomplex_value& _Right)
  {       // construct from float complex value
    _Val[_RE] = _Right._Val[_RE];
    _Val[_IM] = _Right._Val[_IM];
  }

  _Myt& operator=(const _Ty& _Right)
  {       // assign real
    _Val[_RE] = _Right;
    _Val[_IM] = 0;
    return (*this);
  }

  _Myt& operator+=(const _Ty& _Right)
  {       // add real
    _Val[_RE] = _Val[_RE] + _Right;
    return (*this);
  }

  _Myt& operator-=(const _Ty& _Right)
  {       // subtract real
    _Val[_RE] = _Val[_RE] - _Right;
    return (*this);
  }

  _Myt& operator*=(const _Ty& _Right)
  {       // multiply by real
    _Val[_RE] = _Val[_RE] * _Right;
    _Val[_IM] = _Val[_IM] * _Right;
    return (*this);
  }

  _Myt& operator/=(const _Ty& _Right)
  {       // divide by real
    _Val[_RE] = _Val[_RE] / _Right;
    _Val[_IM] = _Val[_IM] / _Right;
    return (*this);
  }

  _Myt& operator+=(const _Myt& _Right)
  {       // add complex
    _Val[_RE] = _Val[_RE] + (_Ty)_Right.real();
    _Val[_IM] = _Val[_IM] + (_Ty)_Right.imag();
    return (*this);
  }

  _Myt& operator-=(const _Myt& _Right)
  {       // subtract complex
    _Val[_RE] = _Val[_RE] - (_Ty)_Right.real();
    _Val[_IM] = _Val[_IM] - (_Ty)_Right.imag();
    return (*this);
  }

  _Myt& operator*=(const _Myt& _Right)
  {       // multiply by complex
    _Ty _Rightreal = (_Ty)_Right.real();
    _Ty _Rightimag = (_Ty)_Right.imag();

    _Ty _Tmp = _Val[_RE] * _Rightreal - _Val[_IM] * _Rightimag;
    _Val[_IM] = _Val[_RE] * _Rightimag + _Val[_IM] * _Rightreal;
    _Val[_RE] = _Tmp;
    return (*this);
  }

  _Myt& operator/=(const _Myt& _Right)
  {       // divide by complex
    _Ty _Rightreal = (_Ty)_Right.real();
    _Ty _Rightimag = (_Ty)_Right.imag();

    if (_Myctraits::_Isnan(_Rightreal) || _Myctraits::_Isnan(_Rightimag))
      _Val[_RE] = _Myctraits::_Nanv(_Rightreal), _Val[_IM] = _Val[_RE];
    else if ((_Rightimag < 0 ? -_Rightimag : +_Rightimag)
             < (_Rightreal < 0 ? -_Rightreal : +_Rightreal))
    {       // |_Right.imag()| < |_Right.real()|
      _Ty _Wr = _Rightimag / _Rightreal;
      _Ty _Wd = _Rightreal + _Wr * _Rightimag;
      if (_Myctraits::_Isnan(_Wd) || _Wd == 0)
      {       // make both parts NaN
        _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
        _Val[_IM] = _Val[_RE];
      }
      else
      {       // compute representable result
        _Ty _Tmp = (_Val[_RE] + _Val[_IM] * _Wr) / _Wd;
        _Val[_IM] = (_Val[_IM] - _Val[_RE] * _Wr) / _Wd;
        _Val[_RE] = _Tmp;
      }
    }
    else if (_Rightimag == 0)
    {       // make both parts NaN
      _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
      _Val[_IM] = _Val[_RE];
    }
    else
    {       // 0 < |_Right.real()| <= |_Right.imag()|
      _Ty _Wr = _Rightreal / _Rightimag;
      _Ty _Wd = _Rightimag + _Wr * _Rightreal;
      if (_Myctraits::_Isnan(_Wd) || _Wd == 0)
      {       // make both parts NaN
        _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
        _Val[_IM] = _Val[_RE];
      }
      else
      {       // compute representable result
        _Ty _Tmp = (_Val[_RE] * _Wr + _Val[_IM]) / _Wd;
        _Val[_IM] = (_Val[_IM] * _Wr - _Val[_RE]) / _Wd;
        _Val[_RE] = _Tmp;
      }
    }
    return (*this);
  }

  _Ty real(const _Ty& _Right)
  {       // set real component
    return (_Val[_RE] = _Right);
  }

  _Ty real() const
  {       // return real component
    return (_Val[_RE]);
  }

  _Ty imag(const _Ty& _Right)
  {       // set imaginary component
    return (_Val[_IM] = _Right);
  }

  _Ty imag() const
  {       // return imaginary component
    return (_Val[_IM]);
  }
};

// CLASS double_complex
class double_complex
  : public _Dcomplex_value
{       // complex with double components
public:
  typedef double _Ty;
  typedef double_complex _Myt;
  typedef _Ctraits_double _Myctraits;
  typedef _Ty value_type;
  // Used to get xcomplex 'templates' correct.
  typedef _Dcomplex_value _MyBaseT;

  double_complex(const float_complex& _Right)
  {       // construct from float complex components
    _Val[_RE] = (_Ty)_Right.real();
    _Val[_IM] = (_Ty)_Right.imag();
  }

  double_complex(const _Ty& _Realval = 0, const _Ty& _Imagval = 0)
  {       // construct from double complex components
    _Val[_RE] = _Realval;
    _Val[_IM] = _Imagval;
  }

  double_complex(const _Dcomplex_value& _Right)
  {       // construct from double complex value
    _Val[_RE] = _Right._Val[_RE];
    _Val[_IM] = _Right._Val[_IM];
  }

  _Myt& operator=(const _Ty& _Right)
  {       // assign real
    _Val[_RE] = _Right;
    _Val[_IM] = 0;
    return (*this);
  }

  _Myt& operator+=(const _Ty& _Right)
  {       // add real
    _Val[_RE] = _Val[_RE] + _Right;
    return (*this);
  }

  _Myt& operator-=(const _Ty& _Right)
  {       // subtract real
    _Val[_RE] = _Val[_RE] - _Right;
    return (*this);
  }

  _Myt& operator*=(const _Ty& _Right)
  {       // multiply by real
    _Val[_RE] = _Val[_RE] * _Right;
    _Val[_IM] = _Val[_IM] * _Right;
    return (*this);
  }

  _Myt& operator/=(const _Ty& _Right)
  {       // divide by real
    _Val[_RE] = _Val[_RE] / _Right;
    _Val[_IM] = _Val[_IM] / _Right;
    return (*this);
  }

  _Myt& operator+=(const _Myt& _Right)
  {       // add complex
    _Val[_RE] = _Val[_RE] + (_Ty)_Right.real();
    _Val[_IM] = _Val[_IM] + (_Ty)_Right.imag();
    return (*this);
  }

  _Myt& operator-=(const _Myt& _Right)
  {       // subtract complex
    _Val[_RE] = _Val[_RE] - (_Ty)_Right.real();
    _Val[_IM] = _Val[_IM] - (_Ty)_Right.imag();
    return (*this);
  }

  _Myt& operator*=(const _Myt& _Right)
  {       // multiply by complex
    _Ty _Rightreal = (_Ty)_Right.real();
    _Ty _Rightimag = (_Ty)_Right.imag();
    _Ty _Tmp = _Val[_RE] * _Rightreal - _Val[_IM] * _Rightimag;
    _Val[_IM] = _Val[_RE] * _Rightimag + _Val[_IM] * _Rightreal;
    _Val[_RE] = _Tmp;
    return (*this);
  }

  _Myt& operator/=(const _Myt& _Right)
  {       // divide by complex
    _Ty _Rightreal = (_Ty)_Right.real();
    _Ty _Rightimag = (_Ty)_Right.imag();

    if (_Myctraits::_Isnan(_Rightreal) || _Myctraits::_Isnan(_Rightimag))
    {       // make both parts NaN
      _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
      _Val[_IM] = _Val[_RE];
    }
    else if ((_Rightimag < 0 ? -_Rightimag : +_Rightimag)
             < (_Rightreal < 0 ? -_Rightreal : +_Rightreal))
    {       // |_Right.imag()| < |_Right.real()|
      _Ty _Wr = _Rightimag / _Rightreal;
      _Ty _Wd = _Rightreal + _Wr * _Rightimag;
      if (_Myctraits::_Isnan(_Wd) || _Wd == 0)
      {       // make both parts NaN
        _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
        _Val[_IM] = _Val[_RE];
      }
      else
      {       // compute representable result
        _Ty _Tmp = (_Val[_RE] + _Val[_IM] * _Wr) / _Wd;
        _Val[_IM] = (_Val[_IM] - _Val[_RE] * _Wr) / _Wd;
        _Val[_RE] = _Tmp;
      }
    }
    else if (_Rightimag == 0)
    {       // make both parts NaN
      _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
      _Val[_IM] = _Val[_RE];
    }
    else
    {       // 0 < |_Right.real()| <= |_Right.imag()|
      _Ty _Wr = _Rightreal / _Rightimag;
      _Ty _Wd = _Rightimag + _Wr * _Rightreal;
      if (_Myctraits::_Isnan(_Wd) || _Wd == 0)
      {       // make both parts NaN
        _Val[_RE] = _Myctraits::_Nanv(_Rightreal);
        _Val[_IM] = _Val[_RE];
      }
      else
      {       // compute representable result
        _Ty _Tmp = (_Val[_RE] * _Wr + _Val[_IM]) / _Wd;
        _Val[_IM] = (_Val[_IM] * _Wr - _Val[_RE]) / _Wd;
        _Val[_RE] = _Tmp;
      }
    }
    return (*this);
  }

  _Ty real(const _Ty& _Right)
  {       // set real component
    return (_Val[_RE] = _Right);
  }

  _Ty real() const
  {       // return real component
    return (_Val[_RE]);
  }

  _Ty imag(const _Ty& _Right)
  {       // set imaginary component
    return (_Val[_IM] = _Right);
  }

  _Ty imag() const
  {       // return imaginary component
    return (_Val[_IM]);
  }
};

inline float_complex::float_complex(const double_complex& _Right)
{       // construct from double complex components
  _Val[_RE] = (_Ty)_Right.real();
  _Val[_IM] = (_Ty)_Right.imag();
}


#define _CMPLX(T)      _CMPLX1(T)
#define _CMPLX1(T)     T##_complex
#define _CTR(T)        _CTR1(T)
#define _CTR1(T)       _Ctraits_##T
#define _TMPLT(T)
#define _Ty    float
#include <xcomplex>    /* define all float_complex functions */

#undef _Ty
#define _Ty    double
#undef _XCOMPLEX_
#include <xcomplex>    /* define all double_complex functions */
#undef _Ty

// FUNCTION operator>>
inline istream& operator>>(istream& _Istr, float_complex& _Right)
{       // extract a float_complex
  typedef float_complex _Myt;
  typedef float _Ty;
  char _Ch = 0;
  _Ty _Real = 0;
  _Ty _Imag = 0;

  if (_Istr >> _Ch && _Ch != '(')
  {       // no leading '(', treat as real only
    _Istr.putback(_Ch);
    _Istr >> _Real;
    _Imag = 0;
  }
  else if (_Istr >> _Real >> _Ch && _Ch != ',')
    if (_Ch == ')')
      _Imag = 0;      // (real)
    else
    {       // no trailing ')' after real, treat as bad field
      _Istr.putback(_Ch);
      _Istr.setstate(ios_base::failbit);
    }
  else if (_Istr >> _Imag >> _Ch && _Ch != ')')
  {       // no imag or trailing ')', treat as bad field
    _Istr.putback(_Ch);
    _Istr.setstate(ios_base::failbit);
  }

  if (!_Istr.fail())
    _Right = _Myt((_Ty)_Real, (_Ty)_Imag);  // store valid result
  return (_Istr);
}

inline istream& operator>>(istream& _Istr, double_complex& _Right)
{       // extract a double_complex
  typedef double_complex _Myt;
  typedef double _Ty;
  char _Ch = 0;
  _Ty _Real = 0;
  _Ty _Imag = 0;

  if (_Istr >> _Ch && _Ch != '(')
  {       // no leading '(', treat as real only
    _Istr.putback(_Ch);
    _Istr >> _Real;
    _Imag = 0;
  }
  else if (_Istr >> _Real >> _Ch && _Ch != ',')
    if (_Ch == ')')
      _Imag = 0;      // (real)
    else
    {       // no trailing ')' after real, treat as bad field
      _Istr.putback(_Ch);
      _Istr.setstate(ios_base::failbit);
    }
  else if (_Istr >> _Imag >> _Ch && _Ch != ')')
  {       // no imag or trailing ')', treat as bad field
    _Istr.putback(_Ch);
    _Istr.setstate(ios_base::failbit);
  }

  if (!_Istr.fail())
    _Right = _Myt((_Ty)_Real, (_Ty)_Imag);  // store valid result
  return (_Istr);
}

// FUNCTION operator<<
inline ostream& operator<<(ostream& _Ostr,
                           const float_complex& _Right)
{       // insert a float_complex
  ostringstream _Sstr;
  _Sstr.flags(_Ostr.flags());
//      _Sstr.imbue(_Ostr.getloc());
  _Sstr.precision(_Ostr.precision());
  _Sstr << '(' << real(_Right) << ',' << imag(_Right) << ')';

  string _Str = _Sstr.str();
  return (_Ostr << _Str.c_str());
}

inline ostream& operator<<(ostream& _Ostr,
                           const double_complex& _Right)
{       // insert a double_complex
  ostringstream _Sstr;
  _Sstr.flags(_Ostr.flags());
//      _Sstr.imbue(_Ostr.getloc());
  _Sstr.precision(_Ostr.precision());
  _Sstr << '(' << real(_Right) << ',' << imag(_Right) << ')';

  string _Str = _Sstr.str();
  return (_Ostr << _Str.c_str());
}
_STD_END
#endif /* _COMPLEX_ */

/*
 * Copyright (c) 1992-2009 by P.J. Plauger.  ALL RIGHTS RESERVED.
 * Consult your license regarding permissions and restrictions.
V5.04:0576 */
